library(ggplot2)
library(plyr)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
library(naniar)
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
library(ggpubr)
## Loading required package: magrittr
##
## Attaching package: 'magrittr'
## The following object is masked from 'package:tidyr':
##
## extract
##
## Attaching package: 'ggpubr'
## The following object is masked from 'package:plyr':
##
## mutate
theme_set(theme_pubclean())
# reading in data
# set working directory
full_data <- read.csv("~/Downloads/bimuno/full_data.csv")
# grouping data to average per cage
full_data_group <- group_by(full_data, cage.code)
full_data_nested <- summarise(full_data_group,
cagecode=first(cage.code),
model=first(model),
drug=first(drug),
Weight_W1=mean(Weight_W1),
Weight_W2=mean(Weight_W2),
Weight_W3=mean(Weight_W3),
Weight_W4=mean(Weight_W4),
weight_euth=mean(Euthanisation_Weight),
OF_Distance_totalcm=mean(OF_Distance_totalcm),
OF_Velocity_meancms=mean(OF_Velocity_meancms),
FST_preswim_strug=mean(FST_preswim_strug),
FST_preswim_swim=mean(FST_preswim_swim),
FST_preswim_imm=mean(FST_preswim_imm),
FST_swim_strug=mean(FST_swim_strug),
FST_swim_swim=mean(FST_swim_swim),
FST_swim_imm=mean(FST_swim_imm),
EPM_close_freq=mean(EPM_close_freq),
EPM_close_secs=mean(EPM_close_secs),
EPM_open_freq=mean(EPM_open_freq),
EPM_open_secs=mean(EPM_open_secs),
EPM_open_latency=mean(EPM_open_latency),
EPM_fullopen_freq=mean(EPM_fullopen_freq),
EPM_fullopen_secs=mean(EPM_fullopen_secs),
EPM_fullopen_latency=mean(EPM_fullopen_latency)
)
# remove NA columns
full_data_nested <- full_data_nested[rowSums(is.na(full_data_nested)) != ncol(full_data_nested), ]
# Variable names:
# [1] "Cage.Code" "animal.code" "model" "drug"
# "Weight_W1" "Weight_W2" "Weight_W3" "Weight_W4" "Euthanisation_Weight" "OF_Distance_totalcm" "OF_Velocity_meancms"
# [7] "FST_preswim_strug" "FST_preswim_swim" "FST_preswim_imm" "FST_swim_strug" "FST_swim_swim" "FST_swim_imm"
# [13] "EPM_close_freq" "EPM_close_secs" "EPM_open_freq" "EPM_open_secs" "EPM_open_latency" "EPM_fullopen_freq"
# [19] "EPM_fullopen_secs" "EPM_fullopen_latency"
###This experiment has two Factors (i.e., types of manipulations: Phenotype and Invasiveness), and the experiment has data for all 4 possible combinations of these two Factors. We should try to analyze it as a 2-way ANOVA.
###Plot the data:
full_data_nested$group <- factor(paste(full_data_nested$model,full_data_nested$drug))
boxplot(FST_swim_imm ~ group, data=full_data_nested,cex.axis = 1.2)
stripchart(FST_swim_imm ~ group, data=full_data_nested,
vertical = TRUE, method = "jitter",
pch = 21, col = "maroon", bg = "bisque",
add = TRUE)
mtext("immobility",2,line=2.5,cex=1.5)
###Model the data
p <- aov(FST_swim_imm ~ group, data=full_data_nested)
###Check the assumptions:
plot(p)
#### R identified case 13, 20 and 20 as outliers and have heterogeneity of residuals.. - all from Gin - Red group.
leveneTest(full_data_nested$FST_swim_imm ~ full_data_nested$model * full_data_nested$drug)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 1.1862 0.3328
## 28
# Levene's Test for Homogeneity of Variance (center = median)
# Df F value Pr(>F)
# group 3 1.1862 0.3328
# 28
# so no significant differences in equal variances - therefore heterogeneity
shapiro.test(p$residuals)
##
## Shapiro-Wilk normality test
##
## data: p$residuals
## W = 0.94601, p-value = 0.1109
# Shapiro-Wilk normality test
#
# data: p$residuals
# W = 0.94601, p-value = 0.1109
# both not significant so no vioaltion of assumptions
hist(resid(p))
# redisuals are somewhat normally distributed
## fst immobility
anova(lm(FST_swim_imm ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: FST_swim_imm
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 10786 10786.1 8.0936 0.008213 **
## drug 1 4453 4453.3 3.3416 0.078221 .
## model:drug 1 361 361.1 0.2710 0.606765
## Residuals 28 37315 1332.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_fst_imm <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(FST_swim_imm),
mean = mean(FST_swim_imm),
sd = sd(FST_swim_imm),
se = sd / sqrt(N)
)
# p1 <- ggplot(sum_data_fst_imm, aes(x=model, y=mean, fill=drug)) +
# geom_bar(position=position_dodge(), stat="identity") +
# geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
# width=.2, # Width of the error bars
# position=position_dodge(.9))+
# labs(x = "", fill= "Drug")
# full_data_nested$dist_cat_n[full_data_nested$model == "Gin"] <- 1
# full_data_nested$dist_cat_n[full_data_nested$model == "Whisky"] <- 2
#
# full_data_nested$scat_adj[full_data_nested$drug == "Blue"] <- -0.2
# full_data_nested$scat_adj[full_data_nested$drug == "Red"] <- 0.2
#
# min.mean.sd.max <- function(x) {
# r <- c(min(x), mean(x) - sd(x), mean(x), mean(x) + sd(x), max(x))
# names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
# r
# }
# p1 <- ggplot(full_data_nested, aes(x=model, y=FST_swim_imm, fill=drug))
# p1 + geom_boxplot(
# aes(color = drug), width = 0.5, size = 0.4,
# position = position_dodge(0.8)
# ) +
# geom_dotplot(
# aes(fill = drug, color = drug),
# #trim = FALSE,
# binaxis='y', stackdir='center', dotsize = 0.8,
# position = position_dodge(0.8)
# )+
# scale_fill_manual(values = c("#0000FF", "#FF0000"))+
# scale_color_manual(values = c("#0000FF", "#FF0000"))
#install.packages("Hmisc")
#library(Hmisc)
p_test <- ggplot(full_data_nested, aes(x=model, y=FST_swim_imm, fill=drug)) +
geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_test + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1),
geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
scale_fill_manual(values = c("#0000FF", "#FF0000"))+
ylab("FST Immobility (secs)")+
labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
#
# p1 <- ggplot(full_data_nested, aes(x=model, y=FST_swim_imm, fill=drug)) +
# geom_boxplot(outlier.size=0)+
# labs(x = "", fill= "Drug")+
# geom_jitter(aes(dist_cat_n + scat_adj,FST_swim_imm),
# # position=position_jitter(width=0.1,height=0),
# # alpha=0.6,
# # size=3,
# show.legend = FALSE)+
# ylab("FST Immobility (secs)")+ stat_summary(fun.y=mean, geom="point", size=2, position=0.2) +
# stat_summary(fun.data = mean_se, geom = "errorbar", position=-0.2)
#
# p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
# effect of model - expect to see difference between FSL and FRL in immobility time
# no effect of drug - which is different to the non-averaged data
Results for FST Swimming:
## fst swimming
anova(lm(FST_swim_swim ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: FST_swim_swim
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 2236.1 2236.13 3.9288 0.05736 .
## drug 1 2153.3 2153.32 3.7833 0.06187 .
## model:drug 1 164.3 164.26 0.2886 0.59537
## Residuals 28 15936.7 569.17
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_fst_swim <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(FST_swim_swim),
mean = mean(FST_swim_swim),
sd = sd(FST_swim_swim),
se = sd / sqrt(N)
)
p1<- ggplot(sum_data_fst_swim, aes(x=model, y=mean, fill=drug)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
width=.2, # Width of the error bars
position=position_dodge(.9))+
labs(x = "", fill= "Intervention")
p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
library(Hmisc)
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:dplyr':
##
## src, summarize
## The following objects are masked from 'package:plyr':
##
## is.discrete, summarize
## The following objects are masked from 'package:base':
##
## format.pval, units
p_swim <- ggplot(full_data_nested, aes(x=model, y=FST_swim_swim, fill=drug)) +
geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_swim + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1),
geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
scale_fill_manual(values = c("#0000FF", "#FF0000"))+
ylab("FST Swimming (secs)")+
labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
# no effect
## fst struggling
anova(lm(FST_swim_strug ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: FST_swim_strug
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 22845 22844.5 18.0540 0.0002148 ***
## drug 1 413 413.3 0.3266 0.5722185
## model:drug 1 1012 1012.5 0.8002 0.3786659
## Residuals 28 35430 1265.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_fst_strug <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(FST_swim_strug),
mean = mean(FST_swim_strug),
sd = sd(FST_swim_strug),
se = sd / sqrt(N)
)
p1<-ggplot(sum_data_fst_strug, aes(x=model, y=mean, fill=drug)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
width=.2, # Width of the error bars
position=position_dodge(.9))+
labs(x = "", fill= "Intervention")
p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
p_strug <- ggplot(full_data_nested, aes(x=model, y=FST_swim_strug, fill=drug)) +
geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_strug + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1),
geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
scale_fill_manual(values = c("#0000FF", "#FF0000"))+
ylab("FST Struggling (secs)")+
labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
#model differences - no effect of drug
# attempt to create stacked bar chart of time spent in FST
# Activity <- c("Immobility", "Swimming", "Struggling")
# Model <- c("Gin Blue", "Gin Red", "Whisky Blue", "Whisky Red")
#
# fst_stacked <- data.frame(full_data_nested$model, full_data_nested$drug, full_data_nested$FST_swim_imm, full_data_nested$FST_swim_swim, full_data_nested$FST_swim_strug)
#
# fst_stacked$full_data_nested.model<- as.factor(fst_stacked$full_data_nested.model)
# fst_stacked$full_data_nested.drug<- as.factor(fst_stacked$full_data_nested.drug)
# fst_stacked$full_data_nested.FST_swim_imm<- as.numeric(fst_stacked$full_data_nested.FST_swim_imm)
# fst_stacked$full_data_nested.FST_swim_swim <- as.numeric(fst_stacked$full_data_nested.FST_swim_swim)
# fst_stacked$full_data_nested.FST_swim_strug <- as.numeric(fst_stacked$full_data_nested.FST_swim_strug)
# barplot(fst_stacked, main="Distribution of Amount of Time Spent in FST",
# xlab="Groups", col=c("darkblue","red"),
# legend = rownames(counts))
# library(lattice)
# barchart(fst_stacked, scales = list(x = "free"),
# auto.key = list(title = "Time Spent in FST"), horizontal=FALSE)
#open field total distance
anova(lm(OF_Distance_totalcm ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: OF_Distance_totalcm
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 1786232 1786232 29.8951 7.736e-06 ***
## drug 1 202 202 0.0034 0.95408
## model:drug 1 210284 210284 3.5194 0.07112 .
## Residuals 28 1673000 59750
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_OF_Distance_totalcm <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(OF_Distance_totalcm),
mean = mean(OF_Distance_totalcm),
sd = sd(OF_Distance_totalcm),
se = sd / sqrt(N)
)
p1 <- ggplot(sum_OF_Distance_totalcm, aes(x=model, y=mean, fill=drug)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
width=.2, # Width of the error bars
position=position_dodge(.9))+
labs(x = "", fill= "Intervention")
p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
p_OF <- ggplot(full_data_nested, aes(x=model, y=OF_Distance_totalcm, fill=drug)) +
geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_OF + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1),
geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
scale_fill_manual(values = c("#0000FF", "#FF0000"))+
ylab("OF Total Distance Travelled (cm)")+
labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
# model differences
#open field speed velocity m/s
anova(lm(OF_Velocity_meancms ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: OF_Velocity_meancms
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 22.0110 22.0110 31.0170 5.867e-06 ***
## drug 1 0.0029 0.0029 0.0041 0.94919
## model:drug 1 3.3168 3.3168 4.6739 0.03932 *
## Residuals 28 19.8700 0.7096
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_of_velo <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(OF_Velocity_meancms),
mean = mean(OF_Velocity_meancms),
sd = sd(OF_Velocity_meancms),
se = sd / sqrt(N)
)
p1 <- ggplot(sum_data_of_velo, aes(x=model, y=mean, fill=drug)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
width=.2, # Width of the error bars
position=position_dodge(.9))+
labs(x = "", fill= "Intervention")
p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
# model differences - and model x drug interaction
#EPM full body to enter open arms frequency
anova(lm(EPM_fullopen_freq ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: EPM_fullopen_freq
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 34.031 34.031 8.5221 0.006852 **
## drug 1 0.125 0.125 0.0313 0.860841
## model:drug 1 2.000 2.000 0.5008 0.484981
## Residuals 28 111.812 3.993
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_epm_fullopen_freq <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(EPM_fullopen_freq),
mean = mean(EPM_fullopen_freq),
sd = sd(EPM_fullopen_freq),
se = sd / sqrt(N)
)
p1 <- ggplot(sum_data_epm_fullopen_freq, aes(x=model, y=mean, fill=drug)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
width=.2, # Width of the error bars
position=position_dodge(.9))+
labs(x = "", fill= "Intervention")
p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
p_EPM_full <- ggplot(full_data_nested, aes(x=model, y=EPM_fullopen_freq, fill=drug)) +
geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_EPM_full + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1),
geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
scale_fill_manual(values = c("#0000FF", "#FF0000"))+
ylab("EPM Frequency of Full Entries to Open Arms")+
labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
# only model differences
# EPM closed arms frequency to enter
anova(lm(EPM_close_freq ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: EPM_close_freq
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 118.195 118.195 34.1733 2.773e-06 ***
## drug 1 18.758 18.758 5.4234 0.02731 *
## model:drug 1 2.820 2.820 0.8154 0.37422
## Residuals 28 96.844 3.459
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_epm_closed_freq <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(EPM_close_freq),
mean = mean(EPM_close_freq),
sd = sd(EPM_close_freq),
se = sd / sqrt(N)
)
p1 <- ggplot(sum_data_epm_closed_freq, aes(x=model, y=mean, fill=drug)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
width=.2, # Width of the error bars
position=position_dodge(.9))+
labs(x = "", fill= "Intervention")
p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
#model and drug differences
## make new variable time spent open / time spent open + closed arms
full_data_nested <- mutate(full_data_nested, EPM_proportion = (EPM_open_secs/(EPM_open_secs + EPM_close_secs)))
full_data_nested <- mutate(full_data_nested, EPM_proportion_total = (EPM_open_secs/300))
anova(lm(EPM_proportion ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: EPM_proportion
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 0.01498 0.0149846 0.7037 0.4087
## drug 1 0.01367 0.0136674 0.6418 0.4298
## model:drug 1 0.00773 0.0077304 0.3630 0.5517
## Residuals 28 0.59626 0.0212950
sum_EPM_proportion <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(EPM_proportion),
mean = mean(EPM_proportion),
sd = sd(EPM_proportion),
se = sd / sqrt(N)
)
p1 <- ggplot(sum_EPM_proportion, aes(x=model, y=mean, fill=drug)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
width=.2, # Width of the error bars
position=position_dodge(.9))+
labs(x = "", fill= "Intervention")
p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
p_EPM_prop <- ggplot(full_data_nested, aes(x=model, y=EPM_proportion_total, fill=drug)) +
geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_EPM_prop + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1),
geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
scale_fill_manual(values = c("#0000FF", "#FF0000"))+
ylab("EPM Proportion of Time Spent in Open Arm (%)")+
labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## no sig differences
anova(lm(EPM_proportion_total ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: EPM_proportion_total
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 0.01598 0.0159758 1.1460 0.2935
## drug 1 0.00850 0.0085043 0.6100 0.4413
## model:drug 1 0.00278 0.0027813 0.1995 0.6585
## Residuals 28 0.39033 0.0139405
sum_EPM_proportion_total <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(EPM_proportion_total),
mean = mean(EPM_proportion_total),
sd = sd(EPM_proportion_total),
se = sd / sqrt(N)
)
p1 <- ggplot(sum_EPM_proportion_total, aes(x=model, y=mean, fill=drug)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
width=.2, # Width of the error bars
position=position_dodge(.9))+
labs(x = "", fill= "Intervention")
p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
# EPM open arms frequency to enter
anova(lm(EPM_open_freq ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: EPM_open_freq
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 13.133 13.1328 3.4477 0.07389 .
## drug 1 2.258 2.2578 0.5927 0.44781
## model:drug 1 6.570 6.5703 1.7249 0.19973
## Residuals 28 106.656 3.8092
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_epm_open_freq <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(EPM_open_freq),
mean = mean(EPM_open_freq),
sd = sd(EPM_open_freq),
se = sd / sqrt(N)
)
p1 <- ggplot(sum_data_epm_open_freq, aes(x=model, y=mean, fill=drug)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
width=.2, # Width of the error bars
position=position_dodge(.9))+
labs(x = "", fill= "Intervention")
p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
# no significant differences
# EPM time spent full body in open arms in secs
anova(lm(EPM_fullopen_secs ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: EPM_fullopen_secs
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 1770 1770.1 1.3293 0.2587
## drug 1 1540 1540.1 1.1566 0.2914
## model:drug 1 128 128.0 0.0961 0.7588
## Residuals 28 37286 1331.6
sum_data_epm_fullopen_secs <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(EPM_fullopen_secs),
mean = mean(EPM_fullopen_secs),
sd = sd(EPM_fullopen_secs),
se = sd / sqrt(N)
)
p1 <- ggplot(sum_data_epm_fullopen_secs, aes(x=model, y=mean, fill=drug)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
width=.2, # Width of the error bars
position=position_dodge(.9))+
labs(x = "", fill= "Intervention")
p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
p_EPM_total_open <- ggplot(full_data_nested, aes(x=model, y=EPM_fullopen_secs, fill=drug)) +
geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_EPM_total_open + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1),
geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
scale_fill_manual(values = c("#0000FF", "#FF0000"))+
ylab("EPM Time spent in Open Arms")+
labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
# no significant differences
# EPM time spent in open arms in secs
anova(lm(EPM_open_secs ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: EPM_open_secs
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 1438 1437.82 1.1460 0.2935
## drug 1 765 765.38 0.6100 0.4413
## model:drug 1 250 250.32 0.1995 0.6585
## Residuals 28 35130 1254.64
sum_data_epm_open_secs <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(EPM_open_secs),
mean = mean(EPM_open_secs),
sd = sd(EPM_open_secs),
se = sd / sqrt(N)
)
p1 <- ggplot(sum_data_epm_open_secs, aes(x=model, y=mean, fill=drug)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
width=.2, # Width of the error bars
position=position_dodge(.9))+
labs(x = "", fill= "Intervention")
p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
# no significant differences
# EPM closed arms time spent in seconds
anova(lm(EPM_close_secs ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: EPM_close_secs
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 75 75.03 0.0551 0.8162
## drug 1 800 800.00 0.5871 0.4499
## model:drug 1 1188 1188.28 0.8721 0.3584
## Residuals 28 38152 1362.56
sum_data_epm_closed_secs <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(EPM_close_secs),
mean = mean(EPM_close_secs),
sd = sd(EPM_close_secs),
se = sd / sqrt(N)
)
p1 <- ggplot(sum_data_epm_closed_secs, aes(x=model, y=mean, fill=drug)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
width=.2, # Width of the error bars
position=position_dodge(.9))+
labs(x = "", fill= "Intervention")
p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
# no significant group differences
# EPM latency to enter open arms
full_data_nested$EPM_open_latency <- as.numeric(full_data_nested$EPM_open_latency)
anova(lm(EPM_open_latency ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: EPM_open_latency
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 300.1 300.12 0.3157 0.5786
## drug 1 225.8 225.78 0.2375 0.6298
## model:drug 1 1300.5 1300.50 1.3682 0.2520
## Residuals 28 26614.8 950.53
sum_data_epm_open_latency <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(EPM_open_latency),
mean = mean(EPM_open_latency),
sd = sd(EPM_open_latency),
se = sd / sqrt(N)
)
p1 <- ggplot(sum_data_epm_open_latency, aes(x=model, y=mean, fill=drug)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
width=.2, # Width of the error bars
position=position_dodge(.9))+
labs(x = "", fill= "Intervention")
p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
# no significant differences
# EPM latency to enter open arms with full body
anova(lm(EPM_fullopen_latency ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: EPM_fullopen_latency
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 182.9 182.88 0.2683 0.6085
## drug 1 130.0 130.01 0.1907 0.6656
## model:drug 1 1098.6 1098.63 1.6119 0.2147
## Residuals 28 19084.1 681.57
sum_data_epm_fullopen_latency <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(EPM_fullopen_latency),
mean = mean(EPM_fullopen_latency),
sd = sd(EPM_fullopen_latency),
se = sd / sqrt(N)
)
p1 <- ggplot(sum_data_epm_fullopen_latency, aes(x=model, y=mean, fill=drug)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
width=.2, # Width of the error bars
position=position_dodge(.9))+
labs(x = "", fill= "Intervention")
p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
## no significant differences
anova(lm(Weight_W1 ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: Weight_W1
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 7595.3 7595.3 11.0173 0.002514 **
## drug 1 5671.1 5671.1 8.2262 0.007763 **
## model:drug 1 1785.0 1785.0 2.5893 0.118808
## Residuals 28 19303.1 689.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_Weight_W1 <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(Weight_W1),
mean = mean(Weight_W1),
sd = sd(Weight_W1),
se = sd / sqrt(N)
)
p1 <- ggplot(sum_data_Weight_W1, aes(x=model, y=mean, fill=drug)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
width=.2, # Width of the error bars
position=position_dodge(.9))+
labs(x = "", fill= "Intervention")
p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
p_weight1 <- ggplot(full_data_nested, aes(x=model, y=Weight_W1, fill=drug)) +
geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_weight1 + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1),
geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
scale_fill_manual(values = c("#0000FF", "#FF0000"))+
ylab("Weight Week 1 (g)")+
labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
anova(lm(Weight_W2 ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: Weight_W2
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 18408.0 18408.0 30.8383 6.129e-06 ***
## drug 1 5791.6 5791.6 9.7024 0.00422 **
## model:drug 1 2476.3 2476.3 4.1485 0.05123 .
## Residuals 28 16713.8 596.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_Weight_W2 <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(Weight_W2),
mean = mean(Weight_W2),
sd = sd(Weight_W2),
se = sd / sqrt(N)
)
p1 <- ggplot(sum_data_Weight_W2, aes(x=model, y=mean, fill=drug)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
width=.2, # Width of the error bars
position=position_dodge(.9))+
labs(x = "", fill= "Intervention")
p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
p_weight2 <- ggplot(full_data_nested, aes(x=model, y=Weight_W2, fill=drug)) +
geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_weight2 + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1),
geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
scale_fill_manual(values = c("#0000FF", "#FF0000"))+
ylab("Weight Week 2 (g)")+
labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
anova(lm(Weight_W3 ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: Weight_W3
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 23247.1 23247.1 33.4939 3.248e-06 ***
## drug 1 6146.6 6146.6 8.8560 0.005962 **
## model:drug 1 2269.7 2269.7 3.2701 0.081304 .
## Residuals 28 19433.9 694.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_Weight_W3 <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(Weight_W3),
mean = mean(Weight_W3),
sd = sd(Weight_W3),
se = sd / sqrt(N)
)
p1 <- ggplot(sum_data_Weight_W3, aes(x=model, y=mean, fill=drug)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
width=.2, # Width of the error bars
position=position_dodge(.9))+
labs(x = "", fill= "Intervention")
p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
p_weight3 <- ggplot(full_data_nested, aes(x=model, y=Weight_W3, fill=drug)) +
geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_weight3 + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1),
geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
scale_fill_manual(values = c("#0000FF", "#FF0000"))+
ylab("Weight Week 3 (g)")+
labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
anova(lm(Weight_W4 ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: Weight_W4
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 26363.8 26363.8 41.0895 6.115e-07 ***
## drug 1 6626.9 6626.9 10.3284 0.003288 **
## model:drug 1 2354.7 2354.7 3.6699 0.065666 .
## Residuals 28 17965.3 641.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_Weight_W4 <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(Weight_W4),
mean = mean(Weight_W4),
sd = sd(Weight_W4),
se = sd / sqrt(N)
)
p1 <- ggplot(sum_data_Weight_W4, aes(x=model, y=mean, fill=drug)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
width=.2, # Width of the error bars
position=position_dodge(.9))+
labs(x = "", fill= "Intervention")
p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
p_weight4 <- ggplot(full_data_nested, aes(x=model, y=Weight_W4, fill=drug)) +
geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_weight4 + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1),
geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
scale_fill_manual(values = c("#0000FF", "#FF0000"))+
ylab("Weight Week 4 (g)")+
labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
anova(lm(weight_euth ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: weight_euth
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 25679.4 25679.4 40.0749 7.557e-07 ***
## drug 1 4548.2 4548.2 7.0978 0.01266 *
## model:drug 1 3270.4 3270.4 5.1037 0.03185 *
## Residuals 28 17942.0 640.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
interaction.plot(full_data_nested$model, full_data_nested$drug, full_data_nested$weight_euth)
interaction.plot(full_data_nested$drug, full_data_nested$model, full_data_nested$weight_euth)
sum_data_weight_euth <- ddply(full_data_nested, c("model", "drug"), summarise,
N = length(weight_euth),
mean = mean(weight_euth),
sd = sd(weight_euth),
se = sd / sqrt(N)
)
p1<- ggplot(sum_data_weight_euth, aes(x=model, y=mean, fill=drug)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
width=.2, # Width of the error bars
position=position_dodge(.9))+
labs(x = "", fill= "Intervention")
p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
# p1 <- ggplot(full_data_nested, aes(x=model, y=weight_euth, fill=drug)) +
# geom_boxplot(outlier.size=0)+
# labs(x = "", fill= "Drug")+
# geom_jitter(aes(dist_cat_n + scat_adj,weight_euth),
# # position=position_jitter(width=0.1,height=0),
# # alpha=0.6,
# # size=3,
# show.legend = FALSE)+
# ylab("Weight (g) at Euthanisation")
#
#
# p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))
p_weight_euth <- ggplot(full_data_nested, aes(x=model, y=weight_euth, fill=drug)) +
geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_weight_euth + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1),
geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
scale_fill_manual(values = c("#0000FF", "#FF0000"))+
ylab("Weight (g) at Euthanisation")+
labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
# graph for showing weight increase over the 4 weeks
sum_data_Weight_W1$Week=1
sum_data_Weight_W2$Week=2
sum_data_Weight_W3$Week=3
sum_data_Weight_W4$Week=4
sum_data_weight_euth$Week=5
all_weight_sum <- rbind(sum_data_Weight_W1, sum_data_Weight_W2, sum_data_Weight_W3, sum_data_Weight_W4, sum_data_weight_euth)
library(ggplot2)
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:Hmisc':
##
## subplot
## The following objects are masked from 'package:plyr':
##
## arrange, mutate, rename, summarise
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
w <- ggplot(all_weight_sum, aes(x = Week, y = mean, fill = drug)) +
geom_bar(stat = "identity",position=position_dodge())+
facet_wrap(~model, nrow = 1)+
labs(x = "", fill= "Intervention")
w+scale_fill_manual(values=c("#0000FF", "#FF0000"))
# weight difference
full_data_nested$weight_diff <- full_data_nested$weight_euth - full_data_nested$Weight_W1
anova(lm(weight_diff ~ model * drug, full_data_nested))
## Analysis of Variance Table
##
## Response: weight_diff
## Df Sum Sq Mean Sq F value Pr(>F)
## model 1 5343.2 5343.2 37.5706 1.293e-06 ***
## drug 1 61.9 61.9 0.4351 0.5149
## model:drug 1 223.1 223.1 1.5690 0.2207
## Residuals 28 3982.1 142.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p_weight_diff <- ggplot(full_data_nested, aes(x=model, y=weight_diff, fill=drug)) +
geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_weight_diff + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1),
geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
scale_fill_manual(values = c("#0000FF", "#FF0000"))+
ylab("Weight (g) Gain across the study")+
labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.